Subspace-based blind identification algorithm of ocean underwater acoustic channel for multi-channel FIR filter

ABSTRACT

The disclosure provides a subspace-based blind identification algorithm of an ocean underwater acoustic channel for multi-channel fir filter, which adopts a technical solution that a channel impulse response coefficient is calculated by quadratic minimization. The disclosure has beneficial effects that estimation precision can be met when using a proper number of samples, and especially when a few noise vectors are used for estimating channel parameters, so that calculation amount is greatly reduced.

CROSS REFERENCE TO RELATED APPLICATION

This application claims the priority of Chinese Patent Application No. 202010043424.1, entitled “Subspace-based Blind Identification Algorithm of Ocean Underwater Acoustic Channel for Multi-channel FIR Filter” filed with the Chinese Patent Office on Jan. 15, 2020, which is incorporated herein by reference in its entirety.

TECHNICAL FIELD

The disclosure relates to the technical field of communication, in particular to a subspace-based blind identification algorithm of an ocean underwater acoustic channel for a multi-channel FIR filter.

BACKGROUND

In a traditional equalization algorithm, an equalizer is designed by sending training sequences or applying priori knowledge to achieve equalization. However, for mobile underwater communication equipments, it is often difficult to extract information accurately by using the traditional equalization algorithm due to position change between communicating parties, so that the communicating parties can not communicate normally. Therefore, it is very important to study blind channel identification or blind channel equalization of underwater acoustic channel. Traditional blind equalization algorithms generally use the characteristics of the signal to propose a cost function, and then use various gradient algorithms for blind identification and blind equalization. Channel identification of MIMO (Multi-input Multi-output) system plays an important role in multi-channel communication system, and can be used to eliminate inter-symbol interference (ISI) in digital signals. The channel identification of MIMO system can complete signal recovery with ensuring the independence of multiple sources. At present, there are mainly two kinds of methods for channel identification, the first kind is a method based on Second-Order Statistics (SOS) for realizing the identification of the system based on the second-order moment of observed signals. The implementation of the method is simple, has closed solution at some cases, requires less data, and can be used to processes Gaussian distributed input signals. However, the SOS method has certain requirements on signals and channels, wherein phase attributes of the channels cannot be obtained according to received signals. The other kind is High-Order Statistics (HOS) method, which carries out channel identification by utilizing high-order moment or high-order spectrum of the observed signals. For the HOS method, at most one signal among input signals can be a Gaussian signal, and a large amount of signal data are to be collected to obtain better high-order statistics. A waveform of the underwater acoustic signal can be recovered and a signal-to-noise ratio of the communication signal can be improved by studying the blind channel identification. In recent years, many estimation methods in the research field of the blind equalization are studied as subspace methods by utilizing mutual orthogonality between signal subspace and noise subspace. The methods have an advantage of obtaining analytical solution of the blind identification problem and are a representative kind of blind equalization method based on second-order statistics.

SUMMARY

The disclosure intends to provide a subspace-based blind identification algorithm of an ocean underwater acoustic channel for a multi-channel fir filter, which has beneficial effects that estimation precision can be met when using a proper number of samples, and especially when a few noise vectors are used for estimating channel parameters, so that calculation amount is greatly reduced.

The technical solution adopted by the disclosure comprises:

calculating a channel impulse response coefficient by quadratic minimization according to formula (16):

$\begin{matrix} {{{q(h)} = {{{K{h}^{2}} - {\sum\limits_{i = 0}^{L + K - 1}{{{\overset{\hat{}}{S}}_{i}^{H}H_{K}}}^{2}}} = {{{K{h}^{2}} - {{h^{H}\left( {\sum\limits_{i = 0}^{L + K - 1}{{\overset{\hat{}}{s}}_{j}{\overset{\hat{}}{s}}_{j}^{H}}} \right)}h}} = {{K{h}^{2}} - {h^{H}\overset{˜}{Q}h}}}}},} & (16) \end{matrix}$

wherein ŝ_(i) represents an M(L+1)×(L+K) filter matrix corresponding to an eigenvector

${\overset{\hat{}}{S}}_{i},{{\overset{˜}{Q} = {\sum\limits_{i = 0}^{L + K - 1}{{\overset{\hat{}}{s}}_{i}{\overset{\hat{}}{s}}_{i}^{H}}}};}$

establishing a formula (17) under a constraint of ∥h∥²=1: {tilde over (q)}(h)=h ^(H) {tilde over (Q)}h  (17);

solving a maximum eigenvalue of {tilde over (Q)} which makes the formula (17) maximized;

estimating a correlation matrix R_(x) by replacing an overall average of the correlation matrix of communication signals with a sample estimation due to the overall average of the correlation matrix of communication signals being unknown in actual communication, wherein

${R_{x} = {\frac{1}{KM}{\sum\limits_{i = 0}^{{KM} - 1}{{\overset{\sim}{X}}_{i}{\overset{\sim}{X}}_{i}^{H}}}}},$ and {tilde over (X)}_(i) represents an observed sample value;

eigenvalue decomposition: decomposing the correlation matrix R_(x) into a form of formula (8): R _(x) =Sdiag(λ₀,λ₁, . . . λ_(L+K−1))S ^(H)+σ² GG ^(H)   (8);

wherein, columns of matrix S are expanded into signal subspaces, and columns of matrix G are expanded into noise subspaces; and

a calculation of a quadratic form: solving min[q(h)]=min[h^(T)Qh] according to formula (15), i.e., solving min(Q_(i))=min(ĝ_(i)ĝ_(i) ^(H)), wherein h=Q⁻¹c.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a SIMO system model;

FIG. 2 shows a comparison of average deviation of parameter estimation;

FIG. 3 shows an MSE comparison;

FIG. 4 shows an average deviation of parameter estimation for different noise vector spaces under quadratic constraints;

FIG. 5 shows an average deviation of parameter estimation for different noise vector spaces under linear constraints;

FIG. 6 shows an MSE comparison for different noise vector spaces under quadratic constraints;

FIG. 7 shows an MSE comparison for different noise vector spaces under linear constraints.

DETAILED DESCRIPTION

Hereinafter, the disclosure will be described in detail with reference to specific embodiments.

1. Blind Identification Model of Multi-Channel FIR Filter

It is assumed that discrete signals transmitted by a signal source are s(k), signals observed after s(k) have been transmitted through channels, filtered and processed are defined as x(k), an impulse response of the channel is h(k) including comprehensive characteristics of transmission, reception, processing of channel and the like, h(k) is linear, and a noise n(k) superimposed during transmission is a band-limited stationary process.

The subspace method relies on two assumptions that: (1) the channel impulse response h(k) has a finite support, i.e. the channel can be represented by an FIR filter; (2) multiple measurements are possible within a character transmission period T.

There are several methods to implement multiple measurements, one of which is adopting multiple sensors to sample signals received by the sensors with the period T, another of which is oversampling a signal received by a single sensor, or a comprehensive application of both. Both the method adopting multiple sensors and the oversampling method adopting a single sensor are equivalent, and they are equivalent to a single-input multiple-output system, as shown in FIG. 1.

Input signal sequence {s(k)} is output through different linear channels {h_(m)(k)}, m=1, 2, . . . , M. An output of any channel can be represented as:

$\begin{matrix} {\left. \begin{matrix} {{x_{m}(k)} = {{{h_{m}(k)} \otimes {s(k)}} + {n_{m}(k)}}} \\ {{X(k)} = {{{H(k)} \otimes {S(k)}} + {N(k)}}} \end{matrix} \right\},} & (1) \end{matrix}$

wherein “⊗” represents convolution, H(k)=[h₁(k), h₂(k), . . . , h_(M) (k)]^(T), X(k)=[x₁(k), x₂(k), . . . , x_(M)(k)]^(T), and N(k)=[n₁(k), n₂(k), . . . , n_(M)(k)]^(T).

If the SIMO system is a finite impulse response system (FIR) with a length of L, when an observation length is K, an output of the SIMO system is: X _(K) =H _(K) S _(K) +N _(K)  (2);

wherein each matrix in formula (2) is expressed as:

${X_{k} = \begin{bmatrix} {X(k)} \\ \vdots \\ {X\left( {k + K - 1} \right)} \end{bmatrix}},{H_{K} = \begin{bmatrix} {H\left( {L - 1} \right)} & \ldots & {H(0)} & \ldots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \ldots & {H\left( {L - 1} \right)} & \ldots & {H(0)} \end{bmatrix}},{S_{k} = \begin{bmatrix} {S\left( {k - K + 1} \right)} \\ \vdots \\ {S(k + K + L)} \end{bmatrix}},{{{{and}\mspace{14mu} N_{k}} = \begin{bmatrix} {N(k)} \\ \vdots \\ {N\left( {k + K - 1} \right)} \end{bmatrix}};}$

wherein X_(k) and N_(k) are both K×M observation signal matrices, H_(K) is an M×K×(K+L) filter matrix, and S_(k) is a (K+L)×1 vector.

It has been proved in Document [71] (i.e. [71] Li Y, Ding Z. Blind Channel Identification Based on Second-order Cyclosttationary statistics. ICASSP, 1993, 4:81-84.) that channel parameters of the SIMO system can be identified by using the second-order statistics output by the system, when M channels meet two following conditions:

-   -   (1) x_(m) ₁ (0)≠0, h_(m) ₁ , (L−1)≠0, 1≤m₁, m₂≤M; and     -   (2) no common zero point in {H_(m)(z)}₁ ^(M), wherein H_(m)(z)         is a Z-transformation of {h_(m)(k)}. Thus H_(K) in the         formula (2) is full-rank for any K≥L−1.

2. Subspace-Based Blind Identification of Multi-Channel FIR Filter

The blind identification is a process of estimating the channel impulse response coefficients H_(k) only according to the observed data X_(k). The formula (2) can be rewritten as: X _(k) =H _(k) S _(k) +N _(k)  (3).

The blind identification is based on an MK×MK autocorrelation matrix R_(x) according to Tong's method (see Document [72] L. Tong, G. Xu, T. Kailath. Blind identification and equalization based on second-oreder statistics: a time domain approach. IEEE Trans. Inform. Theory, 1994, 40(3): 340-350.), i. e: R _(x) =E(X _(k) X _(k) ^(H))  (4).

Formula (5) can be built when the formula (3) is substituted into the formula (4) due to an assumption of additive noise being independent from transmitted signal sequence: R _(x) =H _(k) R _(s) H _(k) ^(H) +R _(N)  (5),

wherein R_(s)=E(S_(k)S_(k)) and R_(N)=E(N_(k)N_(k) ^(H)) represents an autocorrelation matrix of a transmitted signal vector S_(k) and that of a noise vector N_(k) respectively. The autocorrelation matrix R₅ of the signal source has dimensions of (K+L)×(K+L) and is assumed to be full-rank but unknown. The noise autocorrelation matrix has dimensions of KM×KM and is assumed as R_(N)=∝²I for simplicity, wherein a represents a variance of N_(k), I represents a KM×KM unit matrix.

λ₀≥λ₁≥ . . . λ_(KM−1) is assumed to be eigenvalues of R_(x). A signal part of R_(x), i.e. H_(k)H_(s)H_(k) ^(H), has a rank of K+L due to R_(s) being full-rank, that is:

$\begin{matrix} {\left\{ \begin{matrix} {{\lambda_{i} > \sigma^{2}}\ ,{i = 0},1,\ldots\mspace{14mu},{L + K - 1}} \\ {{\lambda_{i} = \sigma^{2}},{i = {L + K}},\ldots\mspace{14mu},{{KM} - 1}} \end{matrix} \right..} & (6) \end{matrix}$

Unit eigenvectors corresponding to the eigenvalues λ₀≥λ₁≥ . . . λ_(K+M−1) are denoted as S₀, S₁, . . . , S_(K+M−1), and unit eigenvectors corresponding to eigenvalues λ_(L+K), λ_(L+K+1), . . . λ_(MK−1) are denoted as G₀, G₁, . . . , G_(KM−(L+K)−1), and there is:

$\begin{matrix} \left\{ {\begin{matrix} {S = \left\lbrack {S_{0},S_{1},\ldots\mspace{14mu},S_{L + K - 1}} \right\rbrack} \\ {G = \left\lbrack {G_{0},G_{1},\ldots\mspace{14mu},G_{{KM} - L - K - 1}} \right\rbrack} \end{matrix},} \right. & (7) \end{matrix}$

wherein S has dimensions of MK×(L+K) and G has dimensions of MK×(MK−L−K).

The matrix R_(s) can then be represented as: R _(x) =Sdiag(λ₀,λ₁, . . . λ_(L+K−1))S ^(H)+σ² GG ^(H)   (8);

Wherein columns of matrix S are expanded into signal subspaces and columns of matrix G are expanded into noise subspaces. The noise subspace is an orthocomplement of the signal subspace, and the signal subspace is also a linear space expanded from columns of a filter H_(K). It can be known from the orthogonality of the noise and signal subspaces that the columns of the filter matrix H_(K) are orthogonal to any vector within the noise subspaces, which can be represented as: G _(i) ^(H) H _(K) =O,0≤i≤MK−L−K  (9).

Although the formula (9) is a linear equation with unknown coefficients, it provides a simple method for blind channel identification. The impulse response coefficients of channels can be uniquely determined by the noise subspace under a condition of at most one constant gain difference with appropriate conditions in the following theorem 1.

Theorem 1: assuming that K≥L, a matrix H_(K−1) is column full rank (i.e., rank(H_(K−1))=L+K−1), and H′_(K) is a non-zero filter matrix with same dimensions as H_(K), a necessary and sufficient condition that a numerical range of H′_(K) is included within a numerical range of H_(K) is that H′_(K) is proportional to H_(K).

The theorem 1 provides a method for estimating channel parameters, wherein the identification of channels depends on a particular structure of H_(K) under a condition that R_(s) is unknown and full-rank.

3. Subspace-Based Parameter Estimation

Orthogonal condition formula (9) is a linear function of channel coefficient, but in practice only a sample estimation Ĝ_(i) of the noise eigenvectors can be used for parameter estimation, so that formula (9) should be solved only in the least square sense. A quadratic cost function is defined as below:

$\begin{matrix} {{q(h)}\overset{def}{=}{\sum\limits_{i = 1}^{{MK} - L - K - 1}{{{{\overset{\hat{}}{G}}_{i}^{H}H_{K}}}^{2}.}}} & \left( {10} \right) \end{matrix}$

It is easily derived from the following lemma 1 that q(h) depends on h rather than H_(K). To this end, the following symbols are introduced, and the structures of h and H_(K) are the same therewith.

It is assumed that v⁽⁰⁾, v⁽¹⁾, . . . , v^((M−1)) are M arbitrary K×1 vectors, denoted as MK×1 vectors of v=[v⁽⁰⁾, v⁽¹⁾, . . . , v^((M−1))]^(T), following symbols, defined as an (L+1)×(L+K) filter matrix V_(L+1) ^((l)), are introduced as:

$\begin{matrix} {{V_{L + I}^{(l)} = \begin{bmatrix} v_{0}^{(l)} & v_{1}^{(l)} & \ldots & v_{K - 1}^{(l)} & 0 & \ldots & 0 \\ 0 & v_{0}^{(l)} & v_{1}^{(l)} & \ldots & v_{K - 1}^{(l)} & 0 & \ldots \\ \vdots & \; & \; & \ldots & \; & \; & \vdots \\ 0 & \ldots & 0 & v_{0}^{(l)} & v_{1}^{(l)} & \ldots & v_{K - 1}^{(l)} \end{bmatrix}};} & (11) \end{matrix}$

wherein V_(L+1) is an M(L+1)×(L+K) matrix formed by stacking M filter matrices V_(L+1) ^((l)), (l=0, 1, . . . M−1), and defined as: V _(L+1)=[V _(L+1) ^((0)T) ,V _(L+1) ^((1)T) , . . . V _(L+1) ^((M−1)T)]^(T)  (12).

Lemma 1: the following structural relationship holds: v ^(T) H _(K) =h ^(T) V _(L+1)  (13).

Obviously, the Lemma 1 holds because H_(K) and V_(L+1) are row full rank and the number of columns of v^(T) and that of h^(T) are equal to the number of rows of H_(K) and V_(L+1). According to the Lemma 1, a quadratic norm of a vector Ĝ_(i) ^(H)H_(K) is defined as: ∥Ĝ _(i) ^(H) H _(K)∥² =Ĝ _(i) ^(H) H _(K) Ĝ _(i) =h ^(H) ĝ _(i) ĝ _(i) ^(H) h  (14);

where ĝ_(i) is an M(L+1)×(L+K) filter matrix corresponding to a vector Ĝ_(i). The formula (10) may be represented as:

$\begin{matrix} {{{q(h)} = {h^{T}Qh}},{Q = {\sum\limits_{i = 0}^{{MK} - L - K - 1}{{\overset{\hat{}}{g}}_{i}{{\overset{\hat{}}{g}}_{i}^{H}.}}}}} & (15) \end{matrix}$

According to the Theorem 1, if the quadratic form is constructed from an actual correlation matrix, then an actual channel impulse response coefficient is unique and therefore q(h)=0. Conversely, if only an estimated covariance matrix is available, a rank of the quadratic form is not equal to M(L+1)−1. Therefore, a nontrivial solution h can be solved out by minimizing the cost function q(h) under appropriate conditions. Typical constraints can be linear constraint and quadratic constraint.

The quadratic constraint minimizes q(h) under a constraint of ∥h∥=1, to solve out a eigenvector corresponding to a minimum eigenvalue of the matrix Q as a solution. The quadratic form has a large calculation amount, but has a good stability.

The linear constraint solves variable h of min[q(h)] under a constraint of c^(H)h=1, wherein c is an M(L+1)×1 vector. The solution at this point is proportional to Q⁻¹c. The linear constraint has a relatively simple calculation, but has a poor stability due to a possibility of Q being singular.

4. Implementation of the Algorithm

The channel impulse response coefficients can be calculated by quadratic minimization. The quadratic form as shown in the formula (15) may be equivalently represented in a form of a series of signal subspace eigenvectors, such as:

$\begin{matrix} {{{q(h)} = {{{K{h}^{2}} - {\sum\limits_{i = 0}^{L + K - 1}{{{\overset{\hat{}}{S}}_{i}^{H}H_{K}}}^{2}}} = {{{K{h}^{2}} - {{h^{H}\left( {\sum\limits_{i = 0}^{L + K - 1}{{\overset{\hat{}}{s}}_{i}{\overset{\hat{}}{s}}_{i}^{H}}} \right)}h}} = {{K{h}^{2}} - {h^{H}\overset{˜}{Q}h}}}}};} & (16) \end{matrix}$

where ŝ_(i) represents an M(L+1)×(L+K) filter matrix corresponding to an eigenvector

${\overset{\hat{}}{S}}_{i},{\overset{˜}{Q} = {\sum\limits_{i = 0}^{L + K - 1}{{\overset{\hat{}}{s}}_{i}{{\overset{\hat{}}{s}}_{i}^{H}.}}}}$ The minimization of the formula (10) under the constraint of ∥h∥²=1 is equivalent to a maximization of the following formula (17): {circumflex over ( )}q(h)=h ^(H) {circumflex over (Q)}h  (17).

The above formula (17) can be maximized by solving a maximum eigenvalue of {circumflex over (Q)}.

The quadratic form with the noise as parameter has a same solution as the quadratic form with the signal as parameter under the constraint of unit norm, however calculation amount of the latter is less than the former, since the former needs to undergo MK−L−K times of summation operations and the latter only needs to undergo L+K times of summation operations.

The theorem 1 points out that the channel impulse response coefficients can be determined by a noise subspace part of the correlation matrix R_(x). The following inference 1 can be derived according to the theorem 1.

Inference 1: parameters of channels can be determined when only a few vectors of the noise subspace are known under appropriate conditions.

The inference 1 is illustrated and proved below. It is proved that channel parameter vector H can be identified by utilizing a vector of the noise subspace for any K>L starting from a simplest case M=2 (i.e. two channels exist).

The inference 1 can be proved according to the steps of: assuming K>L and G≠0 to be an arbitrary vector of the noise space, i.e. G=[G^((0)T),G^((1)T)]^(T) for two channels M=2. Since G_(i) ^(H)H_(K)=O, 0≤i≤MK−L−K, i.e. the noise subspace is orthogonal to the column space of H_(K), it means: G ⁽⁰⁾(z)H ⁽⁰⁾(z)+G ⁽¹⁾(z)H ⁽¹⁾(z)=0  (18);

where G^((i))(z), i=0, 1 is a z-transformation of a K×1 vector G^((i)), i=0, 1. H⁽⁰⁾(z) differs from H⁽¹⁾(z) with only one constant in Formula (18) is proven by a contradiction analytical method as follows:

It is assumed that H′_(K)=[H′^((0)T), H′^((1)T) is a 2(L+1)×1 vector, H′_(K) represents a 2K×(L+K) filter matrix corresponding to H′. A vector G is orthogonal to column of H′_(K), if: G ⁽⁰⁾(z)H′ ⁽⁰⁾(z)+G ⁽¹⁾(z)H′ ⁽¹⁾(z)=0  (19);

where the condition that the formula (18) and the formula (19) cannot hold simultaneously is G⁽⁰⁾(z)H′⁽⁰⁾(z)+G⁽¹⁾(z)H′⁽¹⁾(z)=0. Since H_(K) is assumed to be full-rank, H⁽⁰⁾(z) and H⁽¹⁾(z) are relatively prime. With H⁽⁰⁾(z) being divided by H′⁽¹⁾(z) and H⁽¹⁾(z) being divided by H′⁽¹⁾(z), and since the formula (18) and the formula (19) have the same dimensions, they must be proportional, that is:

$\begin{matrix} {\left\{ \begin{matrix} {{H^{(0)}(z)} = {\alpha\;{{H^{\prime}}^{(0)}(z)}}} \\ {{H^{(1)}(z)} = {\alpha\;{{H^{\prime}}^{(1)}(z)}}} \end{matrix} \right..} & \left( {20} \right) \end{matrix}$

The inference 1 can be further proved in a manner similar to that described above for M≥2.

The above inference 1 shows that the channel parameter vector H may be identified by using one of vectors in the noise subspaces, which means that H can be determined by using one item in the quadratic form.

In fact, L+K linear equations including M(L+1) unknown parameters H will be generated if G₀ is an MK×1 vector in the noise space of matrix R_(x), the vector G₀ is orthogonal to the column space of H_(K), i.e. G₀ ^(H)H_(K)=O or equivalently as H^(H)g₀=0, where g₀ is an M(L+1)×(L+K) filter matrix associated with G₀. There will be no unique solution since the linear system cannot be determined if L+K<M (L+1). p such linear systems can be constructed by selecting p−1 additional vectors in noise subspaces. The system may be determine if the selected p meets p(L+K)≥M(L+1). The implementation and the process of the simulation experiment of the algorithm will be illustrated hereafter.

(1) Implementation of the Algorithm

Due to an overall average of the correlation matrix of communication signals being unknown in actual communication, it is replaced with a sample estimation. A specific implementation includes a step of:

estimating the correlation matrix R_(x), wherein

${R_{x} = {\frac{1}{KM}{\sum\limits_{i = 0}^{{KM} - 1}{{\overset{\sim}{X}}_{i}{\overset{\sim}{X}}_{i}^{H}}}}},$ and {tilde over (X)}_(i) represents an observed sample value, thus R_(x) can meet a Toeplitz structure;

eigenvalue decomposition: decomposing the correlation matrix R_(x) into a form of formula (8); and

a calculation of a quadratic form: solving min[q(h)]=min[h^(r)Qh] according to formula (15), i.e., solving min(Q_(i))=min(ĝ_(i)ĝ_(i) ^(H)), wherein h=Q⁻¹c.

Simulation Experiment:

The algorithm is validated by using Monte Carlo simulation experiments. The transmitted signals adopt a series of independent 4QAM signals, and virtual channels adopt following channels:

H^((0)T) = [(0.049, 0.359), (0.482, −0.569), (−0.556, 0.587), (1.0, 0.0), (−0.171, 0.061)] H^((1)T) = [(0.443, −0.0364), (1.0, 0.0), (0.921, −0.194), (0.189, −0.208), (−0.087, −0.045)] H^((2)T) = [(−0.211, −0.322), (−0.199, 0.918), (1.0, 0.0), (−0.284, −0.524), (0.136, −0.19)] H^((3)T) = [(0.417, 0.030), (1.0, 0.0), (0.873, 0.145), (0.285, 0.309), (−0.049, −0.19)],

where the number of analog channels are M=4, a width of a time window K=10, an intersymbol interference length is L=4, an additional noise is Gaussian white noise, an output signal-to-noise ratio SNR=25 dB, and the numbers of signals in simulation are respectively 250, 500, 1000, 2000, and 5000. 100 times of simulations are performed by using Monte Carlo for each simulation to calculate an average deviation and a mean square error.

(1) Average deviation: the average deviation is defined as:

${b = {\frac{1}{M\left( {L + 1} \right)}{\sum\limits_{i = 1}^{M{({L + 1})}}{{\left\langle H_{i} \right\rangle - H_{i}}}}}},$

wherein

H_(i)

represents an average of the estimated values of H_(i) under Monte Carlo simulation, and H_(i) represents a i^(th) component of an actual parameter vector.

(2) Mean square error: being defined as an average of the mean square errors for each of all unknown parameters;

(3) It can be seen from FIG. 2 that an average deviation of the channel parameter estimation calculated under the quadratic constraint is smaller than an average deviation of the parameter estimation calculated under the linear constraint, especially under a condition of less numbers of signal samples; and it can also be seen from FIG. 3 that results obtained by using MSE is same as those obtained by using average deviation;

(4) The precision of the parameters estimated by using the noise vectors is related to the number of the noise vectors, it can be seen from FIGS. 4-7 that the more the number of the noise vectors is, the higher the estimation precision is whether the linear constraint or the quadratic constraint are used. In addition, it can be observed that the precision of the estimated channel parameters is also increased when the number of signal samples is large, according to the above simulation diagram, but it doesn't mean that the precision of the estimated channel parameters will be obviously increased with an increase of the number of samples, in which the two are not of a linear relationship. It can be seen that the precision of the estimated channel parameters increases slowly and does not change significantly when the number of samples increases from 1500 to 5000, with reference to FIGS. 2-7. It shows that the estimation precision can ensured when a proper number of samples are used, especially the estimation precision can be met when only few of noise vectors are used to estimate the channel parameters, thereby reducing the amount of computation greatly.

(5) So the parameters of channels can be estimated effectively by using several noise vectors, under a condition of an appropriate number of samples, which fully demonstrates correctness and feasibility of the inference 1.

The foregoing is only exemplary embodiments of the disclosure which is not intended to limit the disclosure in any way. Any simple amendments, equivalent variations and modifications of the embodiments according to the technical spirit of the disclosure are within the scope of the technical solution of the disclosure. 

The invention claimed is:
 1. A subspace-based blind identification algorithm of an ocean underwater acoustic channel for multi-channel FIR filter, comprising: calculating a channel impulse response coefficient by quadratic minimization according to a formula (16): $\begin{matrix} {{{q(h)} = {{{K{h}^{2}} - {\sum\limits_{i = 0}^{L + K - 1}{{{\overset{\hat{}}{S}}_{i}^{H}H_{K}}}^{2}}} = {{{K{h}^{2}} - {{h^{H}\left( {\sum\limits_{i = 0}^{L + K - 1}{{\overset{\hat{}}{s}}_{i}{\overset{\hat{}}{s}}_{i}^{H}}} \right)}h}} = {{K{h}^{2}} - {h^{H}\overset{˜}{Q}h}}}}},} & (16) \end{matrix}$ wherein ŝ_(i) represents an M (L+1)×(L+K) filter matrix corresponding to an eigenvector ${\hat{S}}_{i},{{\overset{˜}{Q} = {\sum\limits_{\;^{i = 0}}^{L + K - 1}{{\overset{\hat{}}{s}}_{i}{\overset{\hat{}}{s}}_{i}^{H}}}};}$ establishing a formula (17) under a constraint of ∥h∥²=1: {tilde over (q)}(h)=h ^(H) {tilde over (Q)}h  (17); solving a maximum eigenvalue of {tilde over (Q)} which makes the formula (17) maximized; estimating a correlation matrix R_(x) by replacing an overall average of the correlation matrix of communication signals with a sample estimation due to the overall average of the correlation matrix of communication signals being unknown in actual communication, wherein ${R_{x} = {\frac{1}{KM}{\sum\limits_{i = 0}^{{KM} - 1}{{\overset{˜}{X}}_{i}{\overset{˜}{X}}_{i}^{H}}}}},$  and {tilde over (X)}, represents an observed sample value; eigenvalue decomposition: decomposing the correlation matrix R_(x) into a form of formula (8): R _(x) =Sdiag(λ₀,λ₁, . . . λ_(L+K−1))S ^(H)+σ² GG ^(H)   (8); wherein, columns of matrix S are expanded into signal subspaces, and columns of matrix G are expanded into noise subspaces; and a calculation of a quadratic form: solving min[q(h)]=min[h^(T)Qh] according to formula (15), and solving min(Q_(i))=min(ĝ,ĝ_(i) ^(H)), wherein h=Q⁻¹c. 